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Abatraat  -  In  this  Too— per?  paper,  a  novel  procedure 
for  generating  an  ARMA  spectral  model  of  a  wide  sense 
stationary  time  series  is  developed.  The  parameters 
of  this  model  are  selected  so  that  they  most  closely 
fit  a  set  of  Yule-Welter  equations  which  are  estimated 
from  a  finite  sec  of  time  series'  observations.  This 
ARMA  modeling  method  has  been  found  to  exhibit  a  spec¬ 
tral  estimation  performance  which  is  cypically  superior 
to  such  alternatives  as  Che  maximum  entropy  (AR)  method, 
classical  Fourier  procedures  (MA),and,  the  Box- Jenkins 
method  (ARMA) . 

One  of  the  principal  features  of  this  spectral 
estimation  method  is  the  elegant  algebraic  structure  of 
the  linear  system  of  equations  which  need  be  solved 
when  finding  the  ARMA  model's  parameters.  This  shift- 
invariant  type  structure  gives  rise  to  an  adaptive 
algorithmic  solution  procedure  whose  computational 
efficiency  is  comparable  to  that  achieved  by  recently 
developed  fast  AR  algorithmic  methods.  The  details  of 
Che  adaptive  ARMA  modeling  procedure  will  be  covered  in 
Part  2  of  this  paper.  These  dual  characteristics  of 
excellent  estimation  performance  and  real  time  adaptive 
implementation  mark  this  method  as  being  a  primary 
spectral  estimation  tool. 

I.  INTRODUCTION 

In  many  Interdisciplinary  applications,  it  is  desired 
to  estimate  the  essential  attributes  of  a  generally 
complex  valued  wlde-sense  stationary  time  series 
{x(n)h  Depending  on  the  specific  nature  of  the  time 
series,  this  characterization  is  often  adequately 
revealed  through  knowledge  of  the  time  series'  associa¬ 
ted  autocorrelation  sequence 

r^(n)  -  E{x(n+m)x*(m) }  n-0,  tl,  t  2,...  (1) 

in  which  E  and  *  denote  the  operations  of  expectation 
and  complex  conjugation,  respectively.  On  the  other 
hand,  the  requisite  characterization  may  often  be 
better  made  in  the  frequency  domain  through  the  spect¬ 
ral  density  function 


- 


rx(n)e 


which  is  recognized  as  being  the  Fourier  transform  of 
the  autocorrelation  sequence.  Either  member  of  this 
transform  pair  conveys  che  total  second-ordsr  statis¬ 
tical  information  relative  to  the  underlying  time 
series.  Frequently,  this  second  order  statistical 
characterization  provides  all  the  information  required 
for  a  given  application  (e.g.,  optimal  Wiener  filtering, 
one-step  prediction,  etc.). 

The  classical  spectral  estimation  problem  is  con¬ 
cerned  with  estimating  the  spectral  density  function 
(2)  from  s  finite  sat  of  time  series  observations. 
Without  loss  of  generality,  chess  observations  will  be 
taken  to  be  the  following  N  contiguous  elements 
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x(l),  x(2) . x(N)  ' - -  (3) 

A  variety  of  procedures  Wve  been  proposed  for  using 
these  observations  to  effect  a  spectral  denalty  esti¬ 
mate.  Invariably,  the  resultant  estimate  will  take  on 
a  rational  model  form  as  expressed  by 

bn+b.e  .  .  .  +  b  e 
S  (eJ“)  -  -2 - i - 9 - 

x  1  +  a.e-l“+  .  .  .  ♦  a  e'^ 

1  P 

B  (ej“)|2 

-  -3—-  -  1  (4) 

Ve  >1 

in  which  the  ak  and  b^  are  referred  Co  as  che  model's 
autoregressive  and  moving  average  coefficients,  res¬ 
pectively.  We  shall  refer  to  this  particular  rational 
form  as  an  autoregressive-moving  average  (ARMA)  model 
of  order  (P,q).  It  is  well  known  that  any  continuous 
spectral  density  can  be  approximated  arbitrarily 
closely  by  this  rational  form  if  the  order  pair  (p,q) 
is  selected  adequately  large.  Thus,  by  imposing  a 
rational  form  on  the  spectral  model,  we  incur  no  real 
loss  in  spectral  representation. 

The  preponderance  of  research  and  application  inter¬ 
est  has  been  focused  on  two  special  cases  of  the  above 
ARMA  model.  They  are  the  moving  average  (MA)  model  in 
which  all  of  the  a^  coefficients  are  set  to  zero,  and, 
the  autoregressive  (AR)  model  for  which  all  of  the  bk 
coefficients  except  bg  are  set  to  zero.  The  spectral 
density  estimate  arising  from  a  MA  model  is  seen  to 
possess  no  poles,  and  as  such  it  is  frequently  referred 
to  as  an  all-zero  model.  Using  similar  reasoning,  the 
AR  model  is  referred  to  as  an  all-pole  model,  and,  the 
general  ARMA  model  is  referred  to  as  a  pole-zero  model. 

Classical  Fourier  approaches  [1]  and  che  perlodogram 
method  [2]  are  procedures  which  ultimately  provide  a 
MA  spectral  density  model.  Similarly,  the  maximum 
entropy  method  and  linear  predictive  coding  are  tech¬ 
niques  chat  result  in  AR  spectral  density  models. 
Undoubtedly,  the  primary  reasons  for  interest  in  speci¬ 
al  case  MA  and  AR  models  lie  in  the  fact  chat  they: 

(i)  are  amenable  to  a  tractable  analysis,  (11)  typical¬ 
ly  provide  adequate  spectral  estimation  performance, 
and  (ill)  give  rise  to  coefficient  selection  procedures 
which  are  implementabls  by  computationally  efficient 
algorithms. 

Despits  this  predisposition  towards  MA  snd  AR  models, 
a  growing  interest  in  ARMA  models  is  evident  [3]-[9]. 
This  is  in  recognition  of  the  fact  that  the  more  ■' 

general  ARMA  model  usually  provides  superior  spectral 
estimation  performance  while  at  the  same  time  require* 
fever  model  parameters  to  achieve  that  behavior.  It  is 
because  of  these  very  factors  chst  a  number  of  ARMA 
modeling  procedures  have  been  proposed.  These  include 
the  Box-Jenklns  maximum  likelihood  method  (3),  whiten-  < 
lng  filter  approaches  [4],  [5],  and,  more  recently, 
Cadzow's  high  performance  method  (61— [9 ] .  This  latter 
method  has  been  found  to  provide  a  spectral  estimation 
performance  which  typically  excels  that  obtained  from 
its  MA,  AR,  and  ARMA  counterparts.  1 
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In  this  paper,  we  first  characterise  the  modeling  of 
a  pure  ARMA  time  aeries.  An  analytical  procedure  is 
presented  for  determining  the  underlying  and  b^  co¬ 
efficients  in  which  che  time  series'  actual  auto¬ 
correlation  element  values  are  used.  This  idealistic 
situation  then  provides  che  justification  for  Intro¬ 
ducing  che  high  performance  method  in  which  the  ARMA 
model's  coefficients  are  estimated  from  time  series 
observations  and  not  from  autocorrelation  values.  It 
is  shown  chat  the  p  autoregressive  afc  coefficients  are 
obtained  by  solving  a  consistent  system  of  p  linear 
equations.  When  using  this  direct  approach,  the  com¬ 
plete  set  of  time  series  observations  (3)  are  incorpor¬ 
ated  to  effect  a  single  spectral  estimate  in  one 
computational  effort.  This  approach  is  typically 
referred  to  as  "block  processing".  Moreover,  by  using 
the  generalized  Levinson  algorithm  [101— [11] ,  it  is 
possible  to  solve  the  above  mentioned  system  of  linear 
equations  in  a  computationally  efficient  manner. 

In  Part  2  of  this  paper,  a  recursive  procedure  is 
developed  In  which  the  ARMA  model's  coefficients  are 
updated  as  each  new  time  series  observation  becomes 
available.  In  this  "clme-update  processing"  mode,  an 
adaptive  form  of  spectral  estimation  is  thereby 
achieved.  One  of- the  particularly  attractive  features 
of  this  time-updating  mode  is  its  computational 
efficiency.  Specifically,  the  p  autoregressive  co¬ 
efficients  (in  actuality  prediction  errors)  are 
optimally  updated  with  each  new  time  series  observation. 
The  number  of  multiplication  and  addition  computations 
required  in  this  updating  is  of  che  order  p.  Thus, 
the  computational  complexity  of  the  high  performance 
ARMA  method  is  competitive  with  recently  developed 
"fast"  AR  methods,  but,  its  spectral  estimation  per¬ 
formance  la  typically  far  superior.  Hie  time- update 
mode  Is  particularly  attractive  in  those  situations  in 
which  che  time  series  being  characterized  Is  a  long 
ongoing  process  and  one  wishes  to  generate  a  time 
evolving  sequence  of  spectral  estimates  in  a  real  Clan 
setting. 


II.  ARMA  TIME  SERIES:  PERFECT  MODELING 

In  this  section,  the  second-order  statistical  char¬ 
acterization  of  an  ARMA  time  series  will  be  presented. 
This  characterization  will  play  a  central  role  in  che 
high  performance  spectral  estimation  procedure  that  Is 
to  be  developed  in  subsequent  sections.  The  time 
series  (x(n)}  is  said  to  be  an  ARMA  time  series  of 
order  (p,q)  if  it  la  generated  according  to  the  causal 
linear  recursive  relationship 

q  P 

x(n)  -  l  b  w(n-k)  -  \  a-x(n-k)  (5) 

k-0  k-I 

in  which  (w(n))  Is  a  zero  mean  white  noise  excitation 
whose  individual  elements  have  variance  one.  It  is 
readily  shown  that  the  spectral  density  corresponding 
to  the  response  time  series  (x(n)>  is  given  by  ex¬ 
pression  (4).  Thus,  there  Is  seen  to  be  an  equivalen¬ 
ce  between  a  rational  spectral  density  modal  and  che 
response  of  a  causal  linear  system  to  a  white  noise 
excitation. 

We  will  now  direct  our  attention  to  developing  a 
systematic  procedure  for  identifying  the  recursive 
system's  autoregressive  coefficients  (i.e.,  che  ag) 
and  moving  average  coefficients  (i.e.,  the  b^)  from 
the  response  time  series'  autocorrelation  elements. 

It  will  be  beneficial  to  consider  separately  the  tasks 
of  identifying  these  two  different  sets  of  coefficients. 

Autoregressive  Coefficient  Identification 

The  autoregressive  coefficients  can  be  determined 
directly  upon  examining  che  autocorrelation  character¬ 
ization  of  recursive  system  (3).  This  is  achieved  by 


first  multiplying  both  sides  of  this  recursive  ex¬ 
pression  by  x*(n-m)  and  then  taking  the  expected  value. 
This  is  found  to  result  in  the  well  known  Yule-Walker 
equations 

P 

£  a^r^m-k)  -  -rx(m)  for  o  >  q+1  (6) 

k-1 

where  It  is  Important  to  note  that  the  lag  parameter  a 
is  here  restricted  to  exceed  the  numerator  order  para¬ 
meter  q.  As  a  side  note,  the  Yule-Walker  equations 
will  involve  the  moving  average  coefficients  b^  In  a 
nonlinear  manner  for  lags  0  <  m  <  q.  The  characteristic 
equations  of  expression  (6)  provide  a  straightforward 
procedure  for  obtaining  tne  ARMA  model's  a^  auto¬ 
regressive  coefficients.  This  formally  entails  expres¬ 
sing  the  first  "t"  Yule-Walker  equations  (i.e.,  q+1  <  a 
<  q+t)  in  the  following  matrix  format 


rx(q)  rx(q_1> 

••  rx(q-p+l) 

*1 

rx(qtl) 

rx(q+l)  rx(q) 

..  rx(q-p+2) 

*2 

. . 

rx(q+2) 

•  - 

a 

P 

• 

rx(q+t-l)  rx(q+t-2) 

_ 

..  rx(q-p+t) 

. 

rx(q+t) 

la  which  the  Integer  t  is  taken  to  be  equal  to  or 
larger  than  the  model's  denominator  order  (I.e.,  t  t  p). 
This  linear  system  of  equations  may  be  compactly 
expressed  as 


(8) 


where  R?_  is  a  t»p  autocorrelation  matrix,  r?  la  a  t*l 
cp  c 

autocorrelation  vector,  and ,  a-  is  che  ARMA  model's  p*l 
autoregressive  coefficient  vector.  In  this  representation 
Che  subscripts  t  and  p  are  appended  to  designate  the  num¬ 
ber  of  Yule-Walker  equations  being  used,  and,  the  ARMA 
model's  denominator  order,  respectively.  Similarly,  the 
superscript  q  depicts  the  ARMA  model's  nusMrator  order. 

To  obtain  the  ARMA  model's  autoregreselve  co¬ 
efficients,  one  then  simply  solves  the  consistent  sys¬ 
tem  of  linear  equations  (8).  Valuable  insight  relative 
to  rational  spectral  density  modeling  is  provided  upon 
closer  examination  of  the  autocorrelation  matrix's 
(i.e.,  R^)  algebraic  structure.  It  is  convenient  to 

express  this  characterization  in  the  following  theorem. 

Theorem  1:  Let  (rx(k))  designate  the  auto¬ 
correlation  sequence  which  is  associated  with 
an  ARMA  time  series  of  order  (p,q).  The 
corresponding  system  of  t  linear  aquations 
in  m  unknowns  ss  specified  by 


R, 


(9) 


has  a  unique  solution  provided  that  m»p  and 
n»q  for  any  value  of  tap.  Moreover,  che  rank 
of  the  c*m  matrix  R”a  is  given  by  min  (m.p.t) 

provided  chat  n>q  and,  by  min  (m.t)  for  Osn<q . 

A  proof  of  this  theorem  will  not  be  given  here,  since 
chase  results  are  implicitly  docisMnted  in  various 
textbooks  and  papers  dealing  with  time  series.  It  is 
Important  to  note  chat  even  if  one  has  perfect  auto¬ 
correlation  knowledge  of  an  ARMA  time  aeries,  Che 
evaluation  of  the  associated  autoregressive  co¬ 
efficients  entails  a  determination  of  the  order  pair 
(p , q ) .  This  ordering  information  Is  implicitly  con¬ 
tained  in  Che  algebraic  structure  of  the  autocorrela¬ 
tion  matrix  R°  ,  and,  can  be  obtained  by  examining  this 
CB 
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structure  for  various  combinations  of  che  nonnegative 
integers  m  and  n. 


Bq(eJ“)Bci*(e-,“')  -  Ap(eJlu)Cs*  (e-3“)  +  Ap*  (eJu,)Cs (eJ“) 

(16) 


Moving  Average  Coefficient  Determination 

To  determine  the  coefficients  associated  with  the 
ARMA  time  series,  it  will  be  beneficial  to  Introduce 
the  causal  Image  of  the  time  series'  autocorrelation 
sequence  as  defined  by 

t*  (n)  •  rx(n)u(n)  -  ^r^ (0)B (n)  (10) 

In  which  u(n)  and  5(n)  denote  the  standard  unit-step 
and  unit-Kronecker  delta  sequences,  respectively.  The 
autocorrelation  sequence  may  be  recovered  from  its 
causal  Image  by  using  che  complex  conjugate  symmetry 
property  of  autocorrelation  sequences  (l.e.,  r^i-n)  * 
r]c*(n)).  This  reconstruction  rule  takes  the  form 

rx(n)  •  rx+  (n)  +  rx+ (-n)*  dl) 

Upon  taking  the  Fourier  transform  of  relationship  (11), 
we  have  che  required  spectral  density  expression 

Sx(eJ“)  -  Sx+(e>)  + 

-  melS^  (eJ“)  ]  (12) 

where  Sx+(e3<")  denotes  the  Fourier  transform  of  the 
causal  image  sequence  (rx+(n)}  ■ 

In  whac  is  to  follow,  a  parametric  procedure  for 
representing  S  +(eJ<‘1)  (and  therefore  Sx  (eJ“>))  will  be 

given.  This  will  first  necessitate  che  Introduction 
of  the  auxiliary  sequence 

P 

c(n)  •  r  +(n)  +  £  a.r  (n-k)  ,  0<nsmax(q,p)  (13) 

*  k«l  K  * 

In  which  the  causal  autocorrelation  elements  as  gene¬ 
rated  by  relationship  (10)  and  the  autoregressive  co¬ 
efficients  as  obtained  upon  solving  the  system  of 
equations  (8)  are  used.  According  to  the  Yule-Walker 
equations  (6)  and  the  causal  image  definition  (10), 
it  Is  seen  that  this  auxiliary  sequence  is  identically 
zero  outside  the  indexing  range  0<n<max(q,p) .  With 
this  in  mind,  the  Fourier  transform  of  relationship 
(13)  is  next  taken  and  results  in 

C  (e>)  -  l  c(n)e'Jun  (14.) 

n-0 


-  (1  + 


P 

l 

n«l 


ane-Jun]S  +  (.j“) 


•  A 

»> 


(e>)Sx+(e>> 


(14b) 


in  which  s  ■  max(q,p).  Upon  solving  this  relationship 
for  S^feJu)  and  substituting  this  solution  into  ex¬ 
pression  (12),  the  desired  ARMA  spectral  density  is 
obtained 


.  C  (eJ“) 
S(e>)  -  — — - — 

x  yJ"> 


Ca*(eJu) 

Ap*CJ“> 


* 


A 


z 


(eJ“)Cg(eJ“)  ♦  Ao(.-l")c/(.J“) 

A  (e>)  A  *(.->“) 

P  P 


(IS) 


In  order  to  determine  the  ARMA  model's  b^  moving 
average  coefficients,  we  next  use  this  relationship  in 
conjunction  with  expression  (4)  to  obtain 


A  spectral  factorization  of  this  expression  will  then 
yield  the  prerequisite  b^  coefficients  (assuming  a 
minimum  phase  Bg(e-)w)). 

In  a unwary,  the  spectral  density  and  the  associated 
a|(  and  b^  coefficients  which  characterize  the  ARMA 
time  series  of  order  (P.q)  may  be  determined  by  follow¬ 
ing  the  four  step  procedure  as  outlined  in  Table  1. 

To  carry  out  this  model  identification  scheme,  it  is 
seen  chat  knowledge  of  the  order  pair  (p.q)  and  the 
q+p+1  autocorrelation  elements  rx(0),  rx(l) , • • • >rx(q+p) 
need  be  available. 


1.  Solve  relationship  (8)  for  the  p  autoregressive 
afc  coefficients.  This  will  require  setting  t>p. 


2.  Generate  che  auxiliary  sequence  c(n)  and  its 
Fourier  transform  using  expressions  (13)  and 
(14a),  respectively. 


3.  The  desired  spectral  density  is  then  given  by 
expression  (15). 


4.  Perform  a  spectral  factorization  of  the  poly¬ 
nomial  Bq(eJ“)E(J(eJ“)  as  given  by  equation  (16) 
to  obtain  the  minimum  phase  choice  of  the  b^ 
coefficients. 


Table  1:  Generation  of  the  spectral  density  and  the 
ARMA  model  parameters  associated  with  a  given 
sec  of  autocorrelation  values. 

III.  HIGH  PERFORMANCE  METHOD  OF  ARMA  SPECTRAL 
MODELING 

It  is  possible  to  adapt  many  of  the  ideas  of 
Section  II  to  achieve  an  ARMA  spectral  estimate  when 
only  the  time  series  observations  (3)are  available 
(and  not  autocorrelation  values).  We  shall  again  treat 
separately  the  cases  of  autoregressive  and  moving 
average  coefficient  determination. 

Autoregressive  Coefficient  Estimation 

To  implement  che  autoregressive  coefficient  selection 
process  as  represented  by  relationship  (8)  it  will  be 
necessary  to  compute  appropriate  autocorrelation  esti¬ 
mates  from  che  given  sec  of  clme  series’  observations. 
The  high  performance  ARMA  method  effects  these  esti¬ 
mates  in  the  guise  of  a  convenient  matrix  format  which 
lands  itself  to  a  particularly  efficient  computational 
realization  [ 6 ] — ( 9 ] .  In  particular,  the  autocorrela¬ 
tion  macrlx  and  vector  required  in  expression  (8)  are 
estimated  according  to 

Rtqp  -  YfX  (17) 

rtq  -  Y*x  (18) 

where  the  dagger  symbol  *  denotes  the  operation  of  com¬ 
plex  conjugate  transposition.  The  (N-p)xp  Toeplltz 


type 

matrix 

X  is  specified  by 

*x(p) 

x(p-l)  .  . 

•  x  (1 )  ' 

x(p+l) 

x(p)  .  . 

•  x(2) 

X  - 

i 

x(N-l) 

x(N-2) 

x(N-p) 

while  the  (N-p)»t  Toeplltz  type  matrix  Y  has  the  form 
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i 


x(p-q)  x(p-q-l)  .  .  .  x(p-q-t+l)' 
x(p-q+l)  x(p-q)  .  .  .  x(p-q-C+2) 


appropriate  autoregressive  coefficient  vector,  we  shall 
Introduce  the  following  quadratic  functional 


Y  - 


x(N-q-l)  x(N-q-2)  .  .  .  x(N-q-t) 


and  x  is  an(N-p)xl  vector  given  by* 
x-  [x(p+l),x(p+2),  .  .  .  x(N)]' 


(20)  *  «+A  *  (25) 

in  which  A  is  a  c*t  posltive-semideflnite  diagonal  aatrix 
with  diagonal  elements  A^  that  is  introduced  in 
order  to  provide  one  with  the  option  of  weighting 
differently  the  varioua  error  vector  coaponents.  It  is 
a  simple  natter  to  show  that  an  autoregressive  co¬ 
efficient  vector  which  will  render  this  quadratic 

(21)  functional  a  minima  must  satisfy 


In  formulating  matrix  Y,  we  have  uaed  the  convention 
of  setting  to  zero  any  elements  x(k)  for  which  k  lies 
outside  the  observation  index  range  1  S  k  i  N. 

If  the  autocorrelation  matrix  and  vector  estimates 
(17)  and  (18),  respectively,  are  substituted  into  the 
Yule-Ualker  relationship  (8),  however,  it  is  generally 
found  that  the  resultant  system  of  t  equations  in  the 
p  autoregressive  coefficients  is  inconsistent  for  t>p. 
This  is  due  to  inevitable  inaccuracies  in  the  auto¬ 
correlation  estimates,  and,  to  a  possible  improper  ARMA 
model  order  choice.  In  any  case,  the  system  of 
equations  with  these  estimate  substitutions  will  give 
rise  to  the  t*l  Yule-Walker  approximation  error  vector 
as  specified  by 

a  -  Y+X  a  +  Yfx  (22) 


Upon  caking  the  expected  value  of  e,  it  is  found  that 
for  the  ARMA  modeling  order  choice  in  which  q>p,  that 
this  expectation  results  in 


E{e(k) )  -  (N-q-k) 


rx(q+k)  +  l  anrx(q+k-m) 


lSkst 


(23) 


while  for  the  modeling  order  case  q<p  this  expectation 
produces 

r  r  p 


E(e(k) } 


(N-p)  r  (q+k)  +  l  *  r(q+k-m) 

m-1 


(N-q-k) |rx(q+k)  +  \  amrx(q+k-m) 


.  lSk<p-q 

,  p-q<k<t 
(24) 


In  either  ordering  case,  it  is  seen  that  when  the  tins 
series  is  an  ARMA  process  of  order  (p,q),  the  expected 
value  of  the  error  vector  e  can  be  made  equal  to  zero 
by  a  proper  choice  of  the  autoregressive  coefficient 
vector  a.  Namely,  this  selection  would  be  such  that 
the  underlying  Yule-Walker  equations  (6)  are  satisfied. 
This  implies  that  the  system  of  equetions  (22)  with 
e  -  8  provides  an  unbiased  and  a  consistent  estimate  of 
the.  Yule-Ualker  equations  (8) ,  where  6_  is  the  zero  vector. 

With  the  above  thoughts  in  mind,  an  appealing 
approach  to  selecting  the  autoregressive  coefficient 
vector  is  Immediately  suggested.  Namely,  a  is  chosen 
so  as  to  make  the  error  vector  "as  close”  to  its 
expected  value  of  8.  ae  possible.  This  is  of  course 
predicated  on  the  asaumptlon  that  the  time  series  is  an 
ADM  process  of  order  (p,q)  or  less.  In  order  to 
attain  a  tractable  procedure  for  selecting  an 


*A  more  generalized  version  of  this  estimation  scheme 
can  be  obtained  by  substituting  the  Integer  k  for  p 
wherever  p  appears  in  relationship  (19)-(21).  For  ease 
of  presentation,  k  is  here  restricted  to  be  p. 

2a  little  thought  will  convince  oneself  chat  this  same 
conclusion  will  be  reached  if  both  q  and  p  are  at 
least  equal  to  the  maserator  and  denominator  orders, 
respectively,  of  the  underlying  ASHA  time  series. 


X^YAY^X  a*  --XfYAY+x  (26) 

One  Chen  simply  solves  this  consistent  system  of  p 
linear  equations  in  Che  p  unknown  autoregressive  co¬ 
efficients  to  obtain  an  estimate  for  the  denominator 
of  Che  ARMA  model. 

Moving-Average  Coefficient  Estimation 

There  exist  several  procedures  for  estimating  Che 
ARMA  model's  moving  average  coefficients.  We  shall 
now  briefly  describe  two  procedures  which  have  pro¬ 
vided  satisfactory  performance  and  in  a  sense  comple¬ 
ment  one  another. 


(i)  cfc  Method 


The  procedure  which  has  provided  the  best  fre¬ 
quency  resolution  behavior  is  a  direct  adaption  of 
the  Cfc  method  as  described  in  Section  II  (see  ref. 
[8]).  In  particular,  using  the  set  of  auto¬ 
regressive  coefficient  estimates  aa  obtained  from 
expression  (26)  and  a  suitable  set  of  auto¬ 
correlation  estimates  fx(n)  for  n-0,1 . max(q.p), 

one  computes  the  Cfc  coefficients  using  expression 
(13).  These  coefficients  are  then  uaed  to  achieve 
the  desired  ARMA  spectral  estimate  when  incorporated 
into  relationship  (14a)  and  ultimately  relationship 
(15).  Although  providing  an  excellent  frequency 
resolution  behavior,  this  procedure  suffers  the 
drawback  of  not  having  a  guaranteed  nonnegative 
definite  spectral  density  estimate3.  It  is  with  this 
in  mind  that  the  following  procedure  was  evolved. 

(11)  Smoothed  Periodogrsm  Method 

In  Che  smoothed  periodogrsm  approach,  one  first 
computes  Che  so-called  "residual  time-series 
elements  according  to  the  relationship  (see  ref. [9]) 


P 

t(n)  «  x(n)  +  [  at,x(n-k)  for  p<n<N  (27) 

k-1  * 


in  which  the  ak*  autoregressive  coefficients  as  ob¬ 
tained  by  solving  expression  (26)  are  incorporated. 
From  this  relationship,  it  is  apparent  that  the 
following  spectral  density  expression  holds 


,  S(e>) 

S  (eJ“)  -  r  7— 
*  |A*(ej“) 


(28) 


If  Sx(eJ“)  is  to  correspond  to  an  ARMA  spectral  model 
of  order  (p,q),  it  is  clear  that  a  qth  order  MA 
spectral  estimate  for  the  residual  spectral  density 
Sc (e'“)  must  be  obtained  and  then  substituted  into 
relationship  (26).  The  smoothed  periodogrsm  has 
bean  found  to  be  a  useful  cool  for  this  purpose. 

In  the  smoothed  perlodoRram  method,  one  first 
partitions  the  computed  residual  elements  (27)  into 


■*This  shortcoming  may  be  superficially  avoided  by 
caking  the  absolute  value  of  the  spectral  estimate. 
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L  segments  each  of  length  q+1  as  specified  by 

e.  (n)  «  e (n+p+l+kd)  0<n<q  (29) 

0;k<L-l 


where  "d"  Is  a  positive  inceger  which  specifies  the 
time  shift  between  adjacent  segments.  These  indivi¬ 
dual  segments  will  overlap  if  d<q  and  will  perfectly 
partition  the  residual  sequence  when  d“q+l.  In  order 
to  Include  only  computed  elements ,  the  relevanc  para¬ 
meters  must  be  selected  so  that  q+p+l+(L-l)dsN. 

Next  the  periodogram  for  aach  of  thee  L  segments  is 
taken  and  these  are  averaged  to  obtain  the  desired  q**1 
order  smoothed  periodogram,  chat  is 

.  |2l 

-jwn 


.  1  L-l 

.<eJ“>4  l 


k-0 


1 

q+1 


l  w(n)e,  (n)e 
n-0 


(30) 


J 


where  w(n)  is  a  window  sequence  that  is  normally 
selected  to  be  rectangular  (l.e.,  w(n)-l  for  0<niq). 
The  required  ARMA  spectral  model  is  then  obtained  by 
substituting  this  approximation  into  relationship 
(28)  thereby  giving 


S  (e3“) 


Ss(eJ“) 

:Ap"(eJ“)|2 


(31) 


It  is  readily  shown  that  the  smoothed  periodogram 
procedure  results  in  a  desired  nonnegative  qth  order 
MA  speccral  density  estimate.  Unfortunately,  its  fre¬ 
quency  resolution  capability  is  generally  not  of  the 
same  quality  as  that  of  the  cfc  method.*  On  the  ocher 
hand,  the  smoothed  periodogram  method  provides  more 
smoothly  behaved  spectral  estimates  which  contain  fewer 
spurious  effects. 

To  summarize,  the  required  ARMA  spectral  model  la 
obtained  by  following  the  systematic  procedure  out¬ 
lined  in  Table  2.  The  numerator  dynamic  estimation 
procedure  to  be  used  will  of  course  depend  on  Che  par¬ 
ticular  characteristic  being  sought  (e.g.,  frequency 
resolution,  smoothness,  etc.). 


1.  Specify  values  for  the  ARMA  model's  order 

parameter  pair  (P.q).  the  Yule-Walker  equation 
parameter  t,  and,  the  weighting  matrix's 
diagonal  elements 


2.  Using  the  time  series  observations 

x(l),x(2) . x(N),  construct  the  matrices 

X,  Y,  and  vector  x  according  to  relationships 
(19),  (20),  and  (21),  respectively. 

| - 

|  3.  Determine  the  model's  autoregressive  co- 

I  efficients  by  solving  relationship  (26) 


4.  The  numerator's  dynamics  are  obtained  by  using 
either  the  (1)  c^  method,  or,  (li)  the  smooth¬ 
ed  periodogram  method. 


Table  2.  Basic  steps  of  the  standard  high  per¬ 
formance  ARMA  speccral  estimation  method: 

The  Block  Processing  Mode. 

The  improved  spectral  estimation  performance  ob¬ 
tained  in  using  this  high  performance  method  over  con¬ 
temporary  ARMA  techniques  such  as  Che  Box-Jenkins 
method  Is,  to  a  large  extent,  a  consequence  of  select¬ 
ing  the  Integer  t  to  be  larger  chan  the  minimal 
value  p.  With  the  corresponding  larger  set  of  Yule- 
Walker  equations  that  are  thereby  being  approximated, 
it  intuitively  follows  chat  the  model's  autoregressive 


coefficients  will  be  less  sensitive  to  autocorrelation 
estimate  errors  which  are  embodied  in  Y+X  and  Y+x  than 
would  be  the  case  if  t  were  sec  to  p  (as  in  the  Box- 
Jenklns  method).  This  anticipated  improvement  in 
spectral  estimation  behavior  when  using  Che  high  per¬ 
formance  method  has  in  fact  been  realized  on  a  rather 
large  number  of  numerical  examples  [6]— [9] .  As  we  will 
see  in  part  2  this  high  performance  method  also  lends 
Itself  to  a  particular  fast  adaptive  implementation 
mode  when  c*p.  With  the  two  attributes  of  Improved 
spectral  estimation  performance  and  computational 
efficiency,  this  new  procedure  promises  to  be  an  import¬ 
ant  spectral  estimation  cool. 

It  is  of  interest  to  note  that  whan  q“0  and  t»p,  the 
high  performance  ARMA  spectral  estimation  method  re¬ 
duces  to  the  well  known  AR  covariance  method.  Moreover, 
upon  letting  t  exceed  p,  the  resultant  set  of 
expanded  AR  Yule-Walker  equation  approximations  will 
typically  result  in  better  spectral  estimates  than  the 
standard  AR  covariance  method.  To  Che  authors  know¬ 
ledge,  this  approach  has  not  been  used  in  the  various  AR 
spectral  estimation  procedures  developed  to  dace. 


IV.  ORDER  SELECTION 

One  of  the  important  considerations  when  using  the 
high  performance  method  is  that  of  selecting  the  ARMA 
model  order  pair  (p.q).  This  selection  process  can  be 
made  by  utilizing  properties  of  the  ARMA  autocorrelat¬ 
ion  matrix  as  outlined  in  Theorem  1.  In  particular, 
one  examines  the  column  rank  behavior  of  the  auto¬ 
correlation  matrix  estlmace 


(32) 


that  is  being  used  in  the  high  performance  method. 

Upon  setting  q-t*p,  it  follows  that  the  p«p  auto¬ 
correlation  matrix  estimate  Rp  will  start  becoming 
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ill-conditioned  when  the  order  parameter  p  exceeds  the 
time  series'  inherent  order  value  (assuming  that  q<p) . 
Thus,  the  model  order  determination  can  be  achieved  by 
investigating  the  conditioning  of  Che  matrix  Rpp  as  a 

function  of  p.  As  p  is  Increased,  an  appropriate 
choice  will  be  a  value  p  for  which  there  is  a  precipi¬ 
tate  decrease  in  matrix  conditioning  for  p-p  +  1.  This 
approach,  as  applied  to  the  high  performance  method  of 
spectral  estimation,  has  been  used  successfully  by  Pao 
and  Lee  [13]. 

There  exist  many  matrix  conditioning  measures  which 
may  be  used  for  this  order  determination.  One  of  the 
more  effective  measures  is  the  normalized  determinant 


as  specified  by  ^ 

C (A)  -  det(A)  J  l  l  ia  |‘ 

'  i-1  J-l  13 


(33) 


where  det(A)  designates  the  determinant  of  the  p«p 
matrix  A.  It  is  to  be  noted  chat  this  normalized 
determinant  will  be  zero  when  the  rank  of  A  is  less 
than  p. 


V.  THE  DOWN  SHIFT  OPERATOR 

In  the  analysis  to  follow,  extensive  use  of  the  down 
shift  operator  S  is  made.  This  operstor  down¬ 

shifts  by  one  unit  the  elements  of  the  vector  upon 
which  it  operates  and  Inserts  a  zero  into  the  vacated 
first  component  position.  In  other  words,  this  opera¬ 
tion  takes  the  form 

Sx  •  [0,  x(l>,  x(2) . x(N-l)  ]  *  (34a) 

where  the  N'xl  vector  being  operated  upon  is  given  by 


* 

A  similar  approach  shares  the  same  attributes  as  x  •  ( x ( 1) ,  x(2) .  x(N)]  (34b) 

does  the  smoothed  periodogram.  [12] . 
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The  prims  symbol  hare  used  denotes  Che  operation  of 
vector  transposition.  It  is  a  simple  natter  to  show 
that  Che  downshift  operator  has  the  following  NxN 
matrix  representation 

S  -  [e2  :  e3  :  .  .  .  :  eN  :  9]  05) 

in  which  j)  is  Che  N»1  zero  vector  and  e^  designates  Che 
kth  standard  N*1  basis  vector  whose  components  are  all 
zero  except  for  its  ktl:  which  is  one.  If  this  down¬ 
shift  operator  were  applied  sequentially  m  times  to 
the  vector  x,  it  is  clear  that  a  downshift  of  m  units 
results,  that  is 

S“x  -  [0,  0 . 0,  x(l),x(2) . x(N-m)  (36) 


m  zeros 


VI.  PREWINDOW  MODIFICATION 


In  many  spectral  estimation  applications,  it  is 
necessary  Co  update  the  ARMA  model's  coefficients  as 
new  time  series  observations  become  available.  If  this 
is  to  be  achieved  in  real  time,  however,  it  is  general¬ 
ly  not  feasible  to  apply  the  block  processing  imple¬ 
mentation  of  the  high  performance  method  as  outlined  in 
Table  2.  In  Part  2  of  this  paper,  a  computationally 
efficient  algorithm  for  achieving  this  coefficient  up¬ 
dating  is  developed.  In  order  co  facilitate  this  real 
time  recursive  algorithm,  it  is  necessary  co  slightly 
modify  the  constituent  matrices  X  and  Y,  and  the 
vector  x  which  characterize  the  high  performance  mechod. 
These  modifications  provide  the  required  algebraic 
structure  to  render  the  resultant  modified  high  per¬ 
formance  ARMA  modeling  method  amenable  co  a  compu¬ 
tationally  efficient  recursive  solution. 

Although  a  number  of  modifications  are  possible,  we 
shall  only  treat  the  prewindowing  method  in  this 
Section.5  In  the  premodification  mechod,  the  x  vector 
is  modified  to 

£-  [x(l),  x(2) ,  .  .  .,  x(N)  ] '  (37) 


while  the  X  matrix  is  modified  co  the  Nxp  Toeplitz 
type  matrix 


0 

x(l) 

x(2) 


0 

0 

x(l) 

x(2) 


x(N-l)  x(N-2) 


x  (1) 


x(N-p) 


js*  :  s2*  .  .  .  :  s»x] 


(38) 


where  S  is  the  downshift  operator.  Finally,  the  Y 
matrix  la  modified  to  the  Nxt  Toeplitz  type  matrix 


0 

0  ... 

0 

0 

0 

0 

0 

0 

x(l) 

0 

x(2) 

x  (1 ) 

0 

x(l) 

i(N-q-l) 

x(N-q-2) 

x(N-q-t) 

5 The  posewindowing,  and,  pre  &  postwindowing  modifica¬ 
tion  methods  are  described  in  the  Appendix. 


[s«+1x  :  sq+2x 


(39) 


Upon  examination  of  these  expressions  for  the  modified 
matrices  X  and  7,  it  is  seen  that  they  possess  a  very 
simple  shift  type  structure.  It  is  this  very  structure 
which  renders  the  prewlndowed  modification  amenable 
to  a  computationally  efficient  adaptive  solution  algor¬ 
ithm.  Furthermore,  it  is  to  be  noted  that  lower 
triangular  pxp  and  pxt  matrices  have  been  added  to  the 
top  of  the  original  X  and  Y  matrices  to  form  the  modi¬ 
fied  X  and  y  matrices,  respectively.  These  augmenting 
lower  triangular  matrices  are  uniquely  specified  so  as 
to  make  the  modified  matrices  Toeplitz  in  structure 
(i.e.,  the  elements  along  any  diagonal  are  all  equal) 
with  zeros  appearing  in  the  upper  right  portion  of  each 
matrix.  It  is  this  specific  structure  which  makes  an 
efficient  recursive  solution  possible.  This  method  is 
referred  to  as  prewindowing  since  the  implicit  assump¬ 
tion  that  x(n)  >0  for  n<0  is  being  made. 

If  these  modifications  are  incorporated  into  expres¬ 
sion  (26),  a  modified  set  of  p  linear  equations  in  the 
p  autoregressive  coefficient  unknowns  is  obtained,  that 
is 

.y+yAytx  i°  ■  -x+yAy+=  (40) 

This  system  of  equations  represents  the  least-squares 
solution  to  the  following  statistical  approximation  of 
the  first  t  Yule-Walker  equations 

e  -  y+y  a  +  y  x  (41) 

The  effectiveness  of  this  approximation  can  be  evaluated 
by  taking  the  expected  value  of  this  relationship.  When 
the  ARMA  model  order  parameters  are  such  that  q<p,  this 
expectation  is  found  to  give 


E{e(k) ) 


q+k  p 

(N-q-k)  l  a  r  (q+k-ra)  +  [ 

m-0  m-q+k+1 


(N-m)ai>rx(q+k-m), 

l<k<p-q 


P 

(N-q-k)  l  a  r  (q+k-m) 
m-0  "  * 


p-q<k<t 

(42a) 


where  ag  -  1.  This  implies  that  the  Yule-Walker  equat¬ 
ion  estimate  (42)  is  biased  in  nature.  As  the  data 
length  N  Increases,  however,  this  estimate  becomes 
asymptotically  unbiased.  For  the  ordering  case  q>p,  the 
expectation  is  found  co  yield 
P 

E{e(k)}»  (N-q-k)  Tar  (q+k-m)  ,  l<k<t  (42b) 

m-0  "  x 


which  is  unbiased  in  nature.  Thus,  the  set  of  linear 
equation  estimates (41;  generally  provides  a  satisfactory 
estimate  for  the  associated  Yule-Walker  equations. 

In  order  to  achieve  the  recursive  update  capability 
as  mentioned  previously,  it  will  be  necessary  to 
"restrict"  the  parameter  t  to  be  p.  Thia  in  turn 
results  in  T+y  balng  a  pxp  matrix.  Whan  this  matrix 
is  invertible,  there  always  exists  s  unique  auto¬ 
regressive  vector  which  will  render  the  arror  vector  co 
be  zero,  that  is 

y+y  a*  -  -r^r  (43) 

The  update  algorithm  to  be  presented  in  Part  2,  In 
effect,  allows  us  to  recursively  obtain  Che  solution  for 
the  N+l  data  length  casa  from  the  solution  to  tha  N 
data  length  cssa  [14].  Unfortunately,  the  restriction 
of  t-p  also  generally  results  In  an  associated 
decrease  In  spectral  estimation  performance  (relative 
to  t>p).  Thus,  in  obtaining  a  computationally  efficient 
update  recursive  algorithm,  an  accompanying  dacrease  in 
spectral  estimation  performance  is  the  price  being  paid. 
One  must  therefore  carefully  consider  the  ramifications 
of  this  tradeoff  in  any  given  application.  It  is  note¬ 
worthy,  however,  that  this  performance  degradation 
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diminishes  as  Che  number  of  time  series  observations  N 
grows. 

VII.  GENERALIZED  LEVINSON  ALGORITHM 

In  Che  high  performance  ARMA  modeling  procedures 
presented  in  Sections  III  and  VI,  the  model's  p  auto¬ 
regressive  coefficients  were  obcained  by  solving  a  sys¬ 
tem  of  p  linear  equations.  In  the  special  case  in 
which  t-p  and  the  p*p  matrix  Y+X  is  nonsingular,  this 
relevant  system  of  equations  (26)  simplifies  to 

Y+X  a°  -  -Y+x  (44) 

where  Che  entries  of  the  matrices  X  and  Y  and  the 
vector  x  are  dependent  on  Che  particular  form  being 
used  (l.e.,  unmodified,  previndowed,  postwlndowed ,  etc.) 

If  standard  matrix  inversion  techniques  such  as  the 
Cholesky  decomposition  method  are  used,  on  the  order  of 
p3  multiplications  and  additions  are  required  to  com¬ 
pute  the  solution  to  relationship  (44).  These  standard 
techniques  are  therefore  3aid  to  possess  a  computation¬ 
al  complexity  of  0(p3).  For  relatively  large  values 
of  p,  this  can  result  in  an  undesirable  computational 
burden.  On  the  other  hand,  if  the  pxp  macrix  Y+X  has 
a  near  Toeplitz  structure,  it  is  possible  to  utilize 
the  generalized  Levinson  algorithm  to  obtain  the 
required  solution  using  far  fewer  computations  (10], 
(11).  Since  the  matrix  Y+X  is  being  used  to  approxi¬ 
mate  the  Toeplitz  matrix  Rd  ,  there  is  good  reason  to 
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anticipate  that  Y  X  might  possess  this  structural 
feature. 

To  measure  the  degree  to  which  Y+X  is  Toeplitz  in 
structure,  it  is  necessary  to  introduce  che  concept  of 
displacement  rank.  The  displacement  rank  a (A)  of  the 


pxp  matrix  A  is  formally  given  by 

a (A)  -  mln(a_(A),  a+(A) ] 

(45a) 

where 

a_(A)  «  rank (A  -  SAs' ] 

(45b) 

a+(A)  «  rank(A  -  SAS] 

(45c) 

in  which  S  is  the  aforementioned  down  shift 
(35).  When  the  matrix  A  is  Toeplitz,  it  is 

operacor 

readily 

shown  that  its  displacement  rank  is  two  (or  less). 

Thus,  a  matrix  whose  displacement  rank  is  near  two  la 
said  Co  be  close  to  Toeplitz  in  structure  and  therefore 
amenable  to  efficient  Inversion  using  the  generalized 
Levinson  algorithm. 

If  the  displacement  rank  of  Che  p*p  macrix  Y'X  is  a, 
it  has  been  shown  that  one  can  use  the  generalized 
Levinson  algorithm  to  solve  expression  (44)  with  a 
corresponding  computational  complexity  of  0(up2)^.  If 
>  is  sufficiently  smaller  than  p,  a  significant  com¬ 
putational  savings  can  be  thereby  realized  relative  to 
standard  matrix  inversion  routines.  Fortunately,  the 
displacement  rank  of  YfX  is  adequately  small  for  che 
unmodified  high  performance  ARMA  modeling  method  and 
its  prewindowed  version  (as  well  as  che  postwlndowed 
and  pre  4  postwlndowed  versions).  This  is  a  direct 
consequence  of  the  fact  that  the  columns  of  matrices 
X  and  Y  are  simply  shifted  versions  on  one  another. 

One  may  readily  show  that  the  displacement  rank  of 
matrix  Y+X  for  each  of  the  high  performance  methods  is 
as  shown  in  Table  3.  Since  these  displacement  ranks 
are  so  small,  it  is  clear  chat  the  generalized  Levinson 
algorithm  may  be  advantageous ly  used  for  solving  che 
linear  syecaa  of  equations  (44). 


6 As  a  byproduct  of  this  solution  procedure,  the  optimal 
autoregressive  coefficient  vectors  for  all  ARMA  models 
of  autoregressive  order  k  are  obtained  for  l<k<p. 


Method 

Displacement 
Rank 
a (Y+X) 

Standard 

u 

Prewindow 

3 

Pos  twlndow 

3 

Pre  4  Postwindow 

2 

Table  3:  Displacement  rank  of  the  matrix  Y+X 

for  the  various  high  performance  ARMA  methods 

When  che  parameter  c  is  allowed  to  Increase  beyond  p 
so  as  to  obtain  an  improved  spectral  estimation  per¬ 
formance,  the  displacement  rank  of  each  of  the  methods 
spelled  out  in  Table  3  increases.  It  is  readily  shown 
that  for  t>p  the  displacement  rank  increases  to 
a2(Y+X)  in  all  cases.  For  example,  the  displacement 
rank  of  the  t*p  matrix  Y+X  for  che  standard  procedure 
increases  to  (4)2-16  and  so  forth.  For  excessively 
large  values  of  p,  it  would  then  be  advantageous  to  use 
the  generalized  Levinson  algorithm  to  solve  relationship 
(44)  when  case  t>p.  The  computational  complexity  there¬ 
by  obtained  would  be  on  the  order  of  a2p2. 

VIII.  NUMERICAL  EXAMPLE 

The  unmodified  ARMA  modeling  method  of  spectral  esti¬ 
mation,  as  presented  in  Section  III,  has  been  found  to 
possess  a  significantly  superior  performance  when  com¬ 
pared  to  3uch  contemporary  alternatives  as  the 
perlodogram,  maximum  entropy,  and,  the  Box-Jenkins 
methods  when  applied  to  "narrow"  band  time  series  (i.e., 
sumsof  sinusoids  in  white  noise  [ 6 ] — [ 9 ]  and  [13]). 

With  this  in  mind,  the  effectiveness  of  both  che  un¬ 
modified  and  modified  ARMA  modeling  procedures  will  now 
be  examined  for  a  "moderately  wide  band"  time  series. 

In  particular,  we  shall  treat  the  time  series  as  recent¬ 
ly  considered  by  Bruzzone  and  Kaveh  (15].  Specifically, 
their  ARMA  time  series  of  order  (4,4)  is  characterized 
by  ,  , 

*k  m  *k  +  0. 5e^  (46a) 

1  2 

where  the  individual  time  series  x^  and  x^  are  generated 
according  to 

\  -  °'4vi  -  °-93xk-z +  €i  (46b) 

4‘-°-54-l  -0-93V2  ♦*£ 

1  2 

in  which  che  c^,  s^,  and  are  uncorrelated  Gaussian 

random  variables  with  zero  mean  and  unit  variance.  It 
then  follows  that  the  spectral  density  characterizing 
time  series  (46)  is  given  by 

S  (<v)  -  |l -0.4e*Jul+0.93«'j2‘Ji 

x  r2 

♦  |l  +0.5e'j“’  +0.93e"j2w|  +0.25  (47) 

Using  the  time  series  description  (46),  twenty 
different  sampled  sequences  each  of  length  64  were 
generated.  These  twenty  observation  sets  were  then  used 
to  test  various  spectral  estimation  methods.  In  Figure 
1,  the  twenty  superimposed  plots  of  the  ARMA  model 
spectral  estimates  of  order  (4,4)  obtained  using  the 
first  iterate  of  the  Box-Jenkins  method,  and,  this 
paper's  unmodified  method  with  Afcjt  *  (0.95)*1")  and 
selections  of  t-4,  8,  and,  20  are  shown.  For  compari¬ 
son  purposes,  che  ideal  spectrum  (47)  is  also  shovn. 

From  these  plots,  two  observations  may  be  made: 

(i)  the  unmodified  method  with  t-4  yields  a  marginally 
better  spectral  estimate  than  the  Box-Jenkins  method, 
and,  (ii)  the  unmodified  spectral  estimates  Improve 
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significantly  as  t  is  Increased  free  the  minimal 
value  4.  This  latter  observation  is  nost  noteworthy 
and  indicates  that  the  incorporation  of  acre  chan  the 
alnlaal  nuaber  of  Yule-Walker  aquations  for  determin- 
lng  Che  A3Q4A  model's  autoregressive  coefficients  has 
Che  anticipated  effect  of  significantly  improving 
spectral  estimation  performance. 

Next,  the  modification  methods  developed  in  Section 
V  and  the  appendix  were  applied  to  these  twenty  differ¬ 
ent  saapled  sequences  of  length  64  to  obtain  ARMA  aodel 
spectral  estimates  of  order  (4,4).  The  resultant 
spectra  are  shown  in  Figure  2  where  it  is  apparent  that 
only  "a  modest"  degradation  in  spectral  estimation  per¬ 
formance  has  accrued  due  to  the  transient  effects  in¬ 
troduced  by  the  modified  methods.  This  is  Indeed 
welcomed  news  given  the  ability  to  Implement  these 
modified  methods  with  exceptionally  fast  algorithms. 

It  is  to  be  noted  that  the  "postaodified”,  and  Che 
"pre  &  poscmodlfled"  methods  are  identical  in  this 
example. 

As  a  final  example,  twenty  different  sampled  sequen¬ 
ces  each  of  length  200  were  generated  according  to 
expression  (46).  With  this  longer  data  length,  it  was 
anticipated  that  an  Improvement  in  spectral  estimation 
performance  would  result.  A  marked  improvement  is  in 
fact  realized  as  is  made  evident  from  Figure  3  where 
Che  ARMA  model  spectral  estimates  of  order  (4,4)  are 
shown  for  the  Box-Jenklns  method  and  the  unmodified 
method  for  selections  of  t  ”4,  8,  and  20. 

IX.  CONCLUSION 

A  computationally  efficient  closed  form  method  of 
ARMA  spectral  estimation  has  been  presented.  It  is 
predicated  on  the  approximation  of  a  sec  of  Yule-Walker 
equation  estimates  which  are  generated  from  a  given  set 
of  time  series  observations.  The  ARMA  model's  auto¬ 
regressive  coefficients  are  determined  by  solving  a 
consistent  system  of  linear  equations.  The  displace¬ 
ment  rank  of  the  matrix  corresponding  to  these 
equations  is  four  thereby  indicating  that  an  efficient 
algorithmic  solution  procedure  is  possible. 

The  spectral  estimation  performance  of  this  ARMA 
modeling  procedure  has  been  empirically  found  to  exceed 
that  of  such  counterparts  as  the  maximum  entropy  and 
Box-Jenkins  methods  (e.g.,  see  refs.  [6]-[9]  &  [13]). 
This  behavior  is  Co  a  large  extent,  a  consequence  of 
Che  fact  that  more  chan  the  minimal  number  of  Yule- 
Walker  equation  estimates  are  being  approximated  to 
obtain  the  resultant  ARMA  model  parameters. 

In  order  to  achieve  an  Improved  computational 
efficiency,  a  prewindowed  modification  of  the  proposed 
ARMA  model  spectral  method  was  next  Introduced.  The 
spectral  estimation  performance  of  this  prewindowed 
version  has  been  found  to  be  of  high  quality  for 
moderate  data  lengths.  As  we  will  see  in  Part  2,  this 
prewindowed  method  may  be  implemented  by  an  adaptive 
update  algorithm  whose  computational  efficiency  is 
comparable  to  chat  achieved  by  recently  developed  LMS 
fast  algorithms. 
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APPENDIX 


I.  Postwindow  Modification 


Following  a  similar  procedure  as  employed  in  Section 
VI,  the  addition  of  an  upper  triangular  matrix  to  the 
lower  portion  of  the  matrices  specified  by  equations 
(19)  aiul  (20)  yields  the  prewindowed  matrices 


x(p)  .  . 

.  .  x(l) 

x(N)  . 

.  .  x(N-p+l) 

0  '  ' 

‘  x(N) 

(Al) 
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x(p-q)  ....  x(p-q-t+l) 
x(p-q+l)  ....  x(p-q-t+2) 

x(N)  . x(N-t+l) 


O  x(N)  •  •  x(N4p-q-C) 


(A2) 


where  X  and  I  are  recognized  as  being  (N*p)  and 
(N*t)  Toeplitz  type  matrices,  respectively.  In  a 
similar  manner,  the  column  vector  x  is  modified  to 

*•  [x(p+l),  . .  x(N),  0, . ■ . ,0]  (A3) 
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p  zeros 


The  displacement  rank  of  Cha  matrix  Y*X  Is  raadlly 
found  to  ba  V.  A  generalized  Levinson  procedura  re¬ 
quiring  a  computational  complexity  of  0  Op  2 )  can  than  ba 
appllad  for  solving  cha  systan  of  aquations 

YfX  a  -  -Yfx  (A4) 

A  mors  computationally  officiant  algorithm  aasoclatad 
with  tha  poatvinduw  modification  has  baan  davalopad  [14]. 
It  la  shown  that  tha  nuabar  of  computations  la  raducad 
to  (p  log  p)  if  p  •  q  whara  p  and  q  ara  cha  denominator 
and  numerator  orders  of  cha  AKMA  modal,  respectively . 


II.  Pre  &  Postwindow  Modification  Method 


tha  combination  of  the  previously  discussed  pre- 
vlndowed  and  poscwindowad  modification  mechods  yields 
tha  pre  &  postwindow  modification  method.  The  matrices 
and  vectors  are  modified  in  cha  following  manner. 


u 

.  .  .  .  u 

x(l> 

'  •  0 

X  -  j  x(p) 

.  .  .•  x(l) 

!  x(n) 

.  .  .  x(N-p) 

!  c 

L 

•  x(N) 

r?  •  • 

!  o  .  . 

. 0 

X  ( 1 ) 

.  ‘  1  0 

i 

y  ■  ' 

i 

x(t)  . 

.  .  .  . ‘ .x(i)  ; 

1  x(N)  . 

t+D  j 

1  Q 

•x(N)  ...x(H+p-q-t)j 

L 

£  -  lx(l),  . 

-■  x(N) ,  0,  . 

J 

•a  0] 

p  zeros 


(A5) 


(A6) 


(A7) 


where  X  and  i  denote  (H+p)*p  and  (N+p)*t  Toeplitz  type 
matrices,  respectively,  and  x  denotes  a  (Nfp)xl  column 

vector. 

It  can  be  shown  that  YfX  is  a  Toeplitz  matrix.  The 
conventional  Levinson  algorithm  may  therefore  be  used 
for  solving  che  Toeplitz  system  of  aquations 

y‘la--/tx  (A8) 

in  which  the  inherent  computational  complexity  is 
0(2p-).7  More  recently,  a  fast  algorithmic  solution 
has  been  developed  which  significantly  reduces  this 
computational  complexity. 
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Figure  1:  Spectral  estimates  of  order  (4,4)  by 

the  Box-Jenkins  method  and  by  the  High 
Performance  method  using  various 
values  for  c.  N  ■  64  data  points  for 
each  estimate. 


The  parameter  c  is  here  taken  to  equal  p. 
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